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Abstract. We present a new numerical dissipation algorithm, which can be 
efficiently used in combination with centered finite-difference methods. We start 
from a formulation of centered finite- volume methods for Numerical Relativity, in 
which third-order space accuracy can be obtained by employing just piecewise- 
linear reconstruction. We obtain a simplified version of the algorithm, which can 
be viewed as a centered finite-difference method plus some 'adaptive dissipation'. 
The performance of this algorithm is confirmed by numerical results obtained 
from 3D black hole simulations. 



1. Introduction 

In a recent paper (Alic et al 2007), we presented a centered finite- volume (CFV) 
method for black-hole simulations in numerical relativity. This method is a variant 
of the well known local-Lax- Friedrichs approach (LLF), which is currently being used 
in computational fluid dynamics (including magneto-hydrodynamics). For a specific 
choice of the parameters, this method can be written as a piecewise-fourth-order finite- 
difference (FD) algorithm plus a piecewise-third-order accurate artificial dissipation, 
with automatically tuned local coefRcient. The piecewise prefix comes from the slope 
limiters that are incorporated in order to deal with shocks or other discontinuities. 

Current black hole simulations in Numerical Relativity use instead centered FD 
algorithms combined with a numerical dissipation term of the Kreiss-Oliger type 
(Gustafson et al 1995). This combination can be interpreted as a single numerical 
scheme with built-in dissipation, which can be tuned by a single parameter. In most 
numerical relativity simulations, where only smooth profiles are dealt with, this has 
shown to be an efficient computational approach. In some black hole simulations, 
however, the required amount of dissipation varies from the inner to the outer regions, 
so this approach is lacking some flexibility. 

Our main point is that, as far as the slope limiters are not required, the FV 
algorithm which we developed can be expressed also as a fourth-order centered FD 
algorithms combined with a local dissipation term which is automatically adapted to 
the requirements of the either interior or exterior black hole regions. 

2. The Centered Finite- Volume Method in a Flux-Splitting Approach 

We consider the Einstein field equations written as a system of balance laws 

dtU + dkF\u)^S{u), (1) 
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where the Fliax terms F and the Source terms S depend algebraically on the array 
of dynamical fields u. We will consider first the one-dimensional case. In a regular 
finite difference grid, we choose the elementary cell to be the interval {Xi_i/2, 
centered in the grid point Xi. The resulting discrete scheme is given by 

= " 1^ [ ^+1/2 " ^-1/2 ] + At 5. . (2) 

This general algorithm requires the prescription of the interface fluxes -Fj±i/2- ^^^^ 
use linear reconstruction: the dynamical fields will be modelled as piecewise linear 
functions in each cell. 

Our CFD method (Alic et al 2007) is based on the flux-splitting approach, in 
which the information is evaluated at the grid nodes, selecting the components of the 
flux that will propagate in each direction. In every grid point, the flux can then be 
splitted into two components 

F^^ = F,± X,u, , (3) 

where A is the maximum characteristic speed at the specific grid point. In this way, 
we will have the freedom to choose a different slope for each component, which will 
allow us to improve the space accuracy. 

Then we will consider at each interface two one-sided (left and right) predictions 
from the neighboring points: 

Ft=Ft + \4^ F^=Fr^,-\a-^„ (4) 

where a is the slope of the selected flux component in the corresponding cell. The 
interface flux is obtained by recombining both predictions 

Fi+i/2 = \{Ft+F^). (5) 



3. Third-order-accurate Dissipation Formula 

Let us write the prescription for the slopes of the flux components generically as 

at =aaf + {l- a) af , ar = 6 af + (1 - b) af , (6) 
where a and h are slope coefficients, and we have noted for short 

^Ft- Ft, , of = F^^ - Ft . (7) 

We determine the specific values of the slope coefficients which allow third order 
accuracy, by inserting the slope formulae in the CFV method. Comparing the final 
algorithm with the standard fourth order finite difference method, one obtains the 
values a = 1/3, 6 = 2/3 for the slope coefficients. If no slope limiters are implemented, 
the derivative of the flux can be expressed in closed form as 

Da,{Fi) = -i— [-F,+2 + 8 F,+i - 8 F,_i + F,_2 ] + Dis{F,) , (8) 
12Aa:; 

where the first part of the formula is just the centered fourth-order FD algorithm and 
the second part is the new dissipation term: 

Dis{Fi) = J' [ Xi+2Ui+2 - 4 Xi+iUi+i + 6 XtUi - 4 Xi^iUi-i + Ai_2Uj-2 ] • (9) 
12A.T 

We test the stability and convergence of the resulting algorithm in the Gauge 
Wave Test (Fig. 1), one of the standard tests for Numerical Relativity (Alcubierre et 
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Figure 1. 3D Gauge- Wave test simulations. Tlie profiles in the left panel are 
plots of the metric component g^x for three resolutions (Aa; = 0.02, 0.01, 0.005) 
after 100 crossing times; the last two almost coincide. The right panel shows the 
local convergence rate, calculated by comparing the two higher resolutions with 
the exact solution, confirming third-order convergence. 



al 2004). The plots show that the resulting amount of dissipation is actually very 
small and confirm third-order convergence. Note, however, that in this case all the 
local A coefficients are equal to one, so the new dissipation term coincides with the 
Kreiss-Oliger one for a specific value of its global coefficient. 

4. The 3D Black Hole 

The dissipation algorithm presented above can be easily extended to the 3D case: 

+ 6 >H,j,k ««J,fe - 4:>'i~l.J,k Ui-l,j.k + K-2,3,k Ut-2,].k ] , (10) 

where is the maximum characteristic speed along the x axis, and analogous formulae 
hold for the the y and z axes. 




Figure 2. Lapse evolution in a 3D black hole simulation (zero shift). The dotted 
line profiles arc plotted every IM. The solid line ones are plotted every 5M, up 
to 35M, before boundary-related features become too important (the boundary is 
just at lOM). 
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Let us consider initial data taken from a Schwarzschild black hole 




(11) 



(isotropic coordinates). We will use the 'stuffed black hole' approach (Arbona etal 
1998), by matching a scalar field interior metric to pT|) (the scalar field will also 
evolve). As gauge conditions we choose a singularity-avoidant slicing of the '1+log' 
type in normal coordinates (zero shift). 

We present in (Fig. 2) a low-resolution simulation (Aa; — O.IM) which proves the 
performance of our numerical method in 3D strong- field scenarios. Even in presence 
of steep gradients, the lapse profiles evolve smoothly. 

The numerical tests shown here have been performed with the Z3 evolution system 
(first order in space and time) written in a flux conservative form (Bona et al 1989). 
The time integration is dealt by the well-known method-of-lines, with a third-order 
Runge-Kutta algorithm. 
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